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Abstract 

We propose a simple model for arch formation in silos. We show that small 
pertubations (such as the thermal expansion of the beads) may lead to giant 
stress fluctuations on the bottom plate of the silo. The relative amplitude 
A of these fluctuations are found to be power-law distributed, as A~ T , r ~ 
1.0. These fluctuations are related to large scale 'static avalanches', which 
correspond to long-range redistributions of stress paths within the silo. 
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Although granular matter is very familiar, it does display extremely interesting and 
unexpected features. For example, the stress distribution below a heap of sand shows a 
counter-intuitive minimum right below the apex of the pile [J], reflecting the presence of 
'arches' (or line of forces) within the pile [@-§J. Another very common geometry is that 
of the silo. In this case, it is known since Janssen in 1895 that the weight W supported 
by the bottom plate of a tall silo is only a small fraction of the total weight of the grains 
contained in the silo - most of the weight is 'absorbed' by the side walls |5|,||. Very recent 
experiments actually reveal an even more striking effect: the 'apparent' weight W depends 
very sensitively on temperature J7|. More precisely, a variation of temperature of a few 
degrees induce rather erratic and abrupt changes of W of the order of a several %, while 
relative change in the size of the grains due to thermal expansion is only 10~ 5 — 10~ 4 ! 
Furthermore, these changes are of both signs and hysteretic: the apparent weight W is not 
a unique function of temperature. Another aspect of this phenomenon is the following: 
repeating the same experiment of filling a silo with precisely the same quantity of beads 
lead to an apparent weight W which fluctuates by more than 20% |7|]. 

It has been long recognized that stress propagation within granular media is strongly 
inhomogeneous, with clear 'stress paths' appearing in birefringence experiments. This leads 
to large stress fluctuations which have been recently investigated experimentally || , numer- 
ically and theoretically fH|, through a simple 'scalar' model of stress propagation, where 
only the zz (vertical) component of the stress tensor is considered. Although oversimpli- 
fied, this model leads to predictions in qualitative and perhaps quantitative agreement with 
observations. From a theoretical point of view, it can be shown that the full tensorial 



stress propagation model can indeed reduce, in some particular limits, to the scalar case 
considered by Liu et al. ||[| . 

The basic idea to explain the 'giant' stress fluctuations induced by temperature changes 
(or any other small perturbation) is that the network of 'stress paths' is extremely sensitive 
to small perturbations. A stress path reaching a side wall (which can be seen as a partly 
absorbing boundary for the stress) might change its direction and rather reach the bottom 



plate, leading to a sudden upward 'jump' in W. Conversely, the apparent weight can jump 
down as a path reaching the bottom plate is transformed into one colliding with the side 
walls. 

The model we use to describe more quantitatively this effect is an extension of the one 
considered in to allow for the formation of arches. For a two dimensional packing, each 
'grain' is labelled by two integer (i, n) giving its horizontal and vertical position. [Extension 
of the model to three dimensions is left for future work]. Each grain supports the weight of 
its two upstairs neighbours, and shares its own load randomly between its two downstairs 
neighbours. The corresponding scalar equation for weight propagation is thus: 

W{i, n) = q + {i-l, n- l)W{i - 1, n - 1) 

+ q-(i + l,n- l)W(i + l,ra- 1) +w (1) 

where w is the weight of the grains, assumed to be constant. q±(i,ri) is the fraction of the 
weight transmitted to the the grain (?± 1, n + 1), and is a random variable between zero and 
one, subject to the mass conservation constraint q+(i,n) + q^.(i,n) = 1. The walls are at 
i = ±L, and are simply defined by the fact that q ± (±L, n) = 0, and q T (±L, n) randomly 
chosen between zero and one; a fraction 1 — q T (±L, n) of the mass is thus 'absorbed' by the 
wall. 

As shown by Liu et al. |§, when q + is uniformly distributed, the resulting asymptotic 
weight distribution is given by P(x = = xexp— x, where (W) is the average weight. 
When q + = or 1, the distribution becomes much broader, decaying only as a power-law 

w-l !. 

We want to include the possibility that when the local shear on a given grain is too 
strong, this grain can 'lift off' from one of its downstairs neighbours and therefore only 
transmit its load to the neighbour which is in the direction of the shear. Since we work with 
a scalar model where shear is a priori absent, we assume that the transmitted shear stress 
from (i, n) to (i ± 1, n + 1) is also proportional to q±(z, n)W(i, n), and thus postulate that if 

q+(i - l,n - l)W{i - l,n - 1) 
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- + 1,71 - l)W{i + l,n - 1) > 7fc c W(i, n) (2) 

(where 7£ c is a certain threshold), then g_(z, n) = = 1 — n), meaning that the link to 
the left of (i, n) is removed. In a symmetric fashion, if 

q-(i + l,n - l)W(i + l,n- 1) 

- g+(z - 1, n - 1) W(i - 1, n - 1) > K c W{i, n) (3) 

then q + (i, n) is set to zero. Note that for 1Z C = 1, the model of Liu et al. is exactly recovered 
(all the links are present). Conversely, for 1Z C = 0, the q± can only take values or 1. 

This rule is interesting because it potentially leads to 'static' avalanches in the following 
sense: if TZ C is slightly decreased, a link somewhere is removed - say g_(z,n). Then it is 
quite probable that the threshold will also be exceeded on site (i+l,n+l), thereby inducing 
q_(i + 1, n + 1) to vanish, and so on. In other words, very long arches can suddenly appear 
or disappear when 1Z C is changed (see Fig. 3 below). 

Physically, 7Z C should grow with the density = na d of the packing (d is the dimension 
of space), where n is the number of grains per unit volume and a the size of the grains. 
The reason is that obviously, contacts are easier to removed in a loose packing. It is thus 
reasonnable to assume that TZ C will grow with temperature, with 

5TZ C 5a 



(4) 



1Z C a 

7Z C presumably also depends on the friction coefficient between the grains, and is expected 
to be reach 1Z C = 1 in the limit of zero friction. 

The total weight on the bottom plate of the silo is obviously defined as: 

L 

W= £ W(i,H) (5) 

i=-L 

where H is the height of the silo. Fig 1 shows W as a function of the threshold 1Z C for a given 



sample (i.e. a given set of q± uniformly chosen between zero and one |Tl[) and averaged over 
many samples (thick line), for a silo of size D = 2L + 1 = 61, H = 610. W is a constant 
for small TZ C , and then rapidly grows with TZ C , but reveals enormous fluctuations around the 



average trend. More precisely, we find that for a small change 57Z C of the threshold, then 
with probability p51Z c <C 1 there is a 'static avalanche' in the sense that W changes by a 
finite amount, which is independent of STZ C , and with probability 1 — p5TZ c , W does not 
vary at all. This picture is very similar to the 'shocks' in Burgers' turbulence |T2| . Thus 
the full statistics of the W{1Z C ) curve can be decomposed into two aspects: the density of 
avalanches, and the amplitude of these avalanches. 

• The density p{lZ c ) of 'avalanches' or 'shocks' was obtained by plotting the probability 
that no shocks occur in an interval of length 5TZ C , which is very well approximated by 
exp — [p5TZ c ] (independent shocks). For q uniformly distributed between zero and one, we 
found that p{7Z c ) is of the form HLf(TZ c ), where / is a function of order 1 with increases 
rather rapidly with 1Z C . Typically, for H = 200 grain size and an aspect ratio of 1, the 
density of shocks is ~ 10000. 

• The relative value of the jump of W when an 'avalanche' occurs, will be denoted as A. 
Interestingly, for an aspect ratio of 1, but independently of TZ C (and thus of the density p), 
A is distributed according to a power-law V (A) cx |A| _T , with r = 1.02 ± 0.03 (see Fig. 2). 
On the largest sample (H = D = 201), this power law is observed on nearly five decades, 
cut-off for A smaller than a certain A* (which decreases with H), and cut-off above A = 1, 
which corresponds to relative change of the weight of 100%. 

From a practical point of view, it means that for H = D = 201 (in grain size units), 
a small change of 51Z C as low as 10 -5 will, with probability p51Z c ~ 10 _1 induce a relative 
change A of the apparent weight A, distributed as A _T up to A ~ 1 ! Since r is small (close 
to 1), huge relative changes of the weight between say 10 -2 — 10 _1 have a probability given 
by log - A « , which is of the order of 0.2. Hence, small perturbations typically induce large 
responses. 

Fig 3 shows the 'stress paths' (i.e. grains on which the load a certain arbitrary value) for 
two very close values thresholds such that W has changed dramatically. A small perturbation 
indeed induces a noticable rearrangment of the stress paths (which could perhaps be directly 



detected by acoustic emission [13[]). This phenomenology is actually very close to the one 



of quantum conductance fluctuations in strongly disordered samples, where small changes 
of the chemical potential or the magnetic field lead to a complete change in the optimally 
conducting paths |f4|| . The log of the conductance then also shows large variations. Another 
similar situation is that of 'directed polymers' in random environments, where small changes 
of the local (disordered) potential, or a small external force, can lead to a large-scale change 



in the ground state conformation [|15H l8|| . In the same family of problems, one should 



actually cite spin-glasses, which have also been argued to be 'chaotic', i.e. very sensitive to 



temperature or disorder changes fl9| . |2C 



An interesting consequence of our model is that in the region where small perturbations 
lead to substantial effects, there is a large fraction of links which take a value q± = 0,1. This 
is turn, according to the arguments of ||, should lead to a power-law distribution W~a for 
the local stress on the bottom plate. Fig 4 shows that the stress distribution for 1Z C = 0.9 
can be fitted to W~ 1,18 . Note however that the arguments of assume no correlations 
between the q± = 0, 1, which is not the case here since broken links appear precisely along 
arches. For 1Z C = 1, the distribution has exponential tails [[J. We thus expect that these 
giant fluctuations tend to disappear when the density of the packing is increased, or when 
the friction between the grains decreases. 

When the aspect ratio ^ increases, the situation changes in the following way. When 
a static avalanche occurs 'far' from the bottom plate, the resulting change of weight is 
screened out by the presence of the (partly) absorbing walls. A simple argument allows one 
to understand that the avalanche size distribution V(A) is now cut-off beyond a certain A max 
which decreases as the aspect ratio increases. This is precisely what we observe numerically. 
Furthermore, one observes a progressive 'localisation' of the stress on a unique path, which, 
for large aspect ratios and small TZ C , carries a finite fraction of W. 

Finally, let us discuss the above phenomenon from a slightly different point of view. 
Instead of fixing the disorder (i.e. the q±) and changing the threshold TZ C , one can fix TZ C 
and change the disorder. This would correspond to repeating the same experiment of filling 
a silo, thereby building different local arrangments of the grains. As expected by a kind 



of ergodicity argument, the resulting distribution of apparent rescaled weight 4^ is quite 
broad. As soon as the aspect ratio exceeds 5, this distribution is very close to a gaussian and 
independent of H, with a relative root mean square of ~ 20 — 30%. This compares rather 
well with the experiments. In the same spirit, hysteretic effects can be accounted for by the 
fact that when a contact is lost and then reformed when the temperature comes back to its 
initial value, the new values of q± have no reason to be the same as before, since the local 
environment has slightly changed. 

We have also investigated this model in the heap (sandpile) geometry. When TZ C = 1, 
the stress distribution at the bottom of the pile is maximum just below the apex of the 
pile. However, when TZ C is less than a certain critical value, arches spanning the size of 
the system suddenly proliferate, and a dip at the center of the pile is created, as observed 
experimentally. The importance of long arches to understand the 'dip' was stressed in 0,|J. 
The present model could be a microscopic foundation of this idea. 

As a conclusion, we have shown that a motivated extension of the model proposed by 
Liu et al. to account for arch formation can explain the rather spectacular giant stress 
fluctuations induced by small temperature changes. The basic idea is that of long-range 
modification of stress paths induced by a small perturbation, which can be seen as 'static 
avalanches'. Actually, our model is rather close to the large class of 'SOC (Self-Organized 
Critical) models [T2T],^], initially proposed to describe true 'dynamical' avalanches in granu- 
lar media. The irony would be that these scale invariant avalanches should rather be looked 
for in static properties ! 

Acknowlegments. We want to thank M.E. Cates, E. Clement, J. Duran, J. Rajchenbach 
and J. Wittmer for many fruitful discussions. 

Figure Captions 

Fig 1: Evolution of the apparent weight W as a function of the threshold for a silo of 
size D — 61, and H = 601. The thick line is an average over 1000 samples, while the thin 
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line is a particular realisation of the disorder. 

Fig 2: Distribution of shock amplitudes A, in log-log coordinates, for 1Z C = 0.6 and 
D = L = 201 (thick line), and D = L = 21 (dotted line). For the larger sample, the 
distribution is very nearly given by A -1 between 10~ 6 and 10 _1 . 

Fig 3: Chaoticity of stress paths. We show the large stress paths for two very close 
values of 7Z C , separated by a shock, and H = L = 201. For the second value of 7Z C , the paths 
have been translated by x — > x + 2 for clarity. As one can see, some large scale features of 
the stress network have changed; the total weight carried by the bottom plate has changed 
W = 0.735 to W = 0.845 ! 

Fig 4: Distribution of the local stress Wi on the bottom plate, for 1Z C = 0.9 and a sample 
of size H = 201, D = 201. The best fit gives an exponent —1.18 ± 0.04 for the exponent. 
This value is not far from the value — | expected from the analysis of Coppersmith et al. 
0, valid when q± are independent and can only take the values or 1. 
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